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Abstract 

> 

We find one-, two-, and three-component solitons of the polar and ferromagnetic (FM) types 

(N ■ 

in the general (non-integrable) model of a spinor (three-component) model of the Bose-Einstein 

o 
\o 
o 



condensate (BEC), based on a system of three nonlinearly coupled Gross-Pitaevskii equations. 
The stability of the solitons is studied by means of direct simulations, and, in a part, analytically, 
using linearized equations for small perturbations. Global stability of the solitons is considered by 
means of the energy comparison. As a result, ground-state and metastable soliton states of the 

•i-H . 

FM and polar types are identified. For the special integrable version of the model, we develop the 
Darboux transformation (DT). As an application of the DT, analytical solutions are obtained that 
display full nonlinear evolution of the modulational instability (MI) of a continuous-wave (CW) 
state seeded by a small spatially periodic perturbation. Additionally, by dint of direct simulations, 
we demonstrate that solitons of both the polar and FM types, found in the integrable system, are 
structurally stable, i.e., they are robust under random changes of the relevant nonlinear coefficient 
in time. 
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I. INTRODUCTION 



The experimental observation of matter-wave solitons of the dark A, bright Q, and gap 
types in Bose-Einstein condensates (BECs) have attracted a great deal of attention to 
;he dynamics of nonlinear matter waves, including such aspects as the soliton propagation 
vortex dynamics [5], interference patterns domain walls in binary BECs 7], and the 



modulational instability (MI) 



One of major developments in BECs was the study of spinor condensates. Spinor BECs 
feature an intrinsic three-component structure, due to the distinction between different hy- 
perfine spin states of the atoms. When spinor BECs are trapped in the magnetic potential, 
the spin degree of freedom is frozen. However, in the condensate held by an optical potential, 
the spin is free, making it possible to observe a rich variety of phenomena, such as spin do- 







mams 101 and textures 



investigated in detail |12| - 



111 ] . Recently, properties of BECs with this degree of freedom were 



22| , experimentally and theoretically. In particular, the MI of a 



continuous-wave (CW) state with constant densities in all the three components (sometimes, 
it is called "driven current"), was explored in Ref. [3]. An important results was reported 



in Ref. 



22|, which demonstrated that, under special constraints imposed on parameters, the 



matrix nonlinear Schrodinger (NLS) equation, which is a model of the one- dimensional (ID) 
spinor BEC in the free space, may be integrable by means of the inverse scattering trans- 
form. For that case, exact single-soliton solutions, as well as solutions describing collisions 
between two solitons, were found. 

In this work, we aim to consider several species of simple one-, two-, and three-component 
soliton solutions, of both the ferromagnetic (FM) and polar types, which can be found in 
an analytical form in the general (non-integrable) version of the spinor model. Dynamical 
stability of the solitons will be studied by means of analytical and numerical methods, and 
global stability of different soliton types will be considered by comparison of their energies, 
which will allow us to identify ground-state and metastable soliton states. For the special 
integrable version of the model, we will present a Darboux transformation (DT) and a family 
of exact solutions generated by it from CW states. If the amplitude of the CW background 
vanishes, these solutions go over into solitons found in Ref. 22]. The family also contains 
new spatially periodic time-dependent solutions, that explicitly describe the full nonlinear 
development of the MI, starting from a CW state seeded with an infinitely small spatially 
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periodic perturbation. In addition, we demonstrate structural stability of the solitons of 
both polar and FM types in the integrable model, by showing that they are robust under 
random changes of a nonlinear coefficient whose value determines the integrability. 

The paper is organized as follows. In Section 2, the model is formulated. Various species 
of exact soliton solutions in the non-integrable model, and their stability, are considered 
in Section 3. Section 4 introduces the DT, presents solutions obtained by means of this 
technique, and also numerical results demonstrating the structural stability of the solitons. 
The paper is concluded by Section 5. 



II. THE MODEL 



We start with the consideration of an effectively one-dimensional BEC trapped in a pencil- 
shaped region, which is elongated in x and tightly confined in the transverse directions y, z 
(see, e.g., Ref. p). Therefore, atoms in the F = 1 hyperfme state can be described by a 
ID vectorial wave function, &(x,t) = t), $o(^> t), ^)] T , with the components 

corresponding to the three values of the vertical spin projection, mp = +1, 0, —1. The wave 
functions obey a system of coupled Gross-Pitaevskii (alias NLS) equations 

h 2 

ihd t $ ±1 = -— <9 2 $ ±1 + (c + c 2 )(|$±i| 2 + |$ | 2 )$±i 
2m 

+ (c - c 2 )|$ T i| 2 $±i + c 2 % 1 ®l, 

h 2 

ihd t ® = -— <9 2 $ + (c + c 2 )(|$+i| 2 + |$_i| 2 )$o 
2m 

+ co|$o| 2< I ) o + 2c 2 $+i$_ 1 $*, (1) 
where Co = (go+2g 2 )/3 and c 2 = (#2— <7o)/3 denote effective constants of the mean-field (spin- 



independent) and spin-exchange interaction, respectively |14|. Here gj = Ah aj /[ma ± (l — 
caf/a±)], with / = 0,2, are effective ID coupling constants, aj is the s-wave scattering 
length in the channel with the total hyperfme spin /, a± is the size of the transverse ground 
state, m is the atomic mass, and c = —£(1/2) ~ 1.46. Redefining the wave function as 
$ — > (4> +1 , "s/2005 0-i) T an d measuring time and length in units of h/\co\ and ^/H 2 /2m\co\, 
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respectively, we cast Eqs. (0) in a normalized form 

zd t ±1 = -<9 2 ±1 - (v + a)(|0 ±1 | 2 + 2|0 O | 2 )0 ±1 
-(z/-a)|0 Tl | 2 ±1 - ^ 2 



|2 



^ + a)(|0 +1 | 2 + |0_!| 2 )0 O - 2a0 +1 0_ 1 0*, (2) 



ff = / rfa;<j|^ + J%|^_J 2 -i(z/ + a)(|0 +1 | 4 + |0„ 1 | 4 )-(z/-a)|0 +1 | 21 -' 2 



-ii 



where i/ = — sgn (co), a = — C2 / 1 Co | . 

These equations can be derived from the Hamiltonian 

-+oo f 2 2 

+2 [\d x <P \ 2 - v |0 O | 4 -(u + a) (|0 +1 | 2 + l^l 2 ) |0 O | 2 - a feCi^o + (<^) 2 )] } > 

(3) 

which is a dynamical invariant of the model (dH/dt = 0). In addition, Eqs. (J2J) conserve 
the momentum, 

/+oo 
(0+1^0+1 + <K-xdx<t>-i + 20o«9x0o) ^ 
-oo 

and the solution's norm, proportional to total number of atoms, 

/+oo 
|0 +1 (a;,t)| 2 + |^ 1 (a;,t)| 2 + 2|0 o (a:,t)| 2 rfx. (4) 
-oo 

An obvious reduction of equations (J2J) can be obtained by setting O = 0. In this case, the 
model reduces to a system of two equations which are tantamount to those describing light 
transmission in bimodal nonlinear optical fibers, with the two modes representing either 
different wavelengths or two orthogonal polarizations |23| : 

id t (j) ±1 = -<9 2 </> ±1 - (u + a)\(j) ±1 \ 2 (j) ±1 - [y - a)\<p Tl \ 2 (p ±l . (5) 

In particular, the case of two linear polarizations in a birefringent fiber corresponds to 
(y — a) I '{y + a) = 2/3, i.e., a = u/5, and two different carrier waves or circular polarizations 
correspond to [y — a) / \v + a) = 2, i.e., a = —u/3. The two nonlinear terms in Eqs. (J5J) 
account for, respectively, the SPM and XPM (self- and cross-phase-modulation) interactions 
of the two waves. The MI of CW (uniform) states in the system of two XPM-coupled 
equations © has been studied in detail |23J]. If the nonlinearity in Eqs. (jHJ) is self-defocusing, 
i.e., v — — 1 and 1 ±a > 0(in other words, \a\ < 1), the single NLS equation would of course 
show no MI, but the XPM-coupled system of NLS equations (0) gives rise to MI, provided 
that the XPM interaction is stronger than SPM, i.e., < a < 1. 



III. EXACT SINGLE-, TWO-, AND THREE-COMPONENT SOLITON SOLU- 
TIONS 



The previous consideration of solitons in the model of the spinor BEC was focused on 
the particular integrable case, a — v — 1 [22[. As follows from the above derivation, this 
special case is physically possible, corresponding to ci = Co, or, equivalently, 2g = — g 2 > 0. 
The latter condition can be satisfied by imposing a condition on the scattering lengths, 

a± = 3ca a 2 /(2a + a 2 ), (6) 

provided that 0002(02 — ao) > holds. 

However, it is also possible to find simple exact soliton solutions of both polar and FM 
types in the general case, a 7^ v. Below, we present several soliton species for this case, and 
study their stability. 



A. Single-component ferromagnetic soliton 

First, a single-component FM soliton is given by a straightforward solution (provided 
that a + v > 0), 

0_! = O = 0, <f> +1 = ,/^5-sech (V=H e"^, (7) 

V v + a v ' 

where the negative chemical potential /1 is the intrinsic parameter of the soliton family. While 
the expression ((7|) corresponds to a zero- velocity soliton, moving ones can be generated from 
it in an obvious way by means of the Galilean transformation. Note that the condition 
a + v > implies a > 1 or a > - lin the cases of the, respectively, repulsive [y = — 1) 
and attractive [y = +1) spin-independent interaction. The norm (0J) and energy (jH)) of this 
soliton are 

N = ^t, H=- {v + af N\ (8) 
v + a ' 48 y J 

The soliton (JJJ) is stable against infinitesimal perturbations, as the linearization of Eqs. 
(J2J about this solution demonstrates that all the three equations for small perturbations are 
decoupled. Then, as the standard soliton solution of the single NLS equation is always stable, 
solution (JZJ) cannot be unstable against small perturbations of <f> +1 . Further, the decoupled 
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linearized equations for small perturbations 0_ 1 and O of the other fields are 

id t 0_! = -d 2 x <\)_ x -{y- a)|0 +1 | 2 0-i, (9) 
= -«9 2 o - (z/ + a)|0 +1 | 2 o , (10) 

and it is well known that such linear equations, with |0 +1 | 2 corresponding to the unperturbed 
soliton (J7|), give rise to no instabilities either. 

B. Single- component polar soliton 

The simplest polar soliton, that has only the O component, can be found for v — +1 
(note that the solution does not depend on the other nonlinear coefficient a): 

0o = v / -y"sech (a/^/Ix) e" iM , <p ±l = 0. (11) 

The stability problem for this soliton is more involved, as the unperturbed field O gives rise 
to a coupled system of linearized equations for infinitesimal perturbations ±1 of the other 
fields: 

(id t + n) x±i = ~dlx±i + 2/isech 2 {yf^x) [(1 + a)x ±1 + ax^J , (12) 
where we have substituted v — +1 and redefined the perturbation, 

±1 = exp (-i/it) x±v (13) 
The last term in Eqs. ()12J) gives rise to a parametric gain, which may be a well-known source 



of instability (see, e.g., Ref. 2^). This instability can be easily understood in qualitative 
terms if one replaces Eqs. (fT2|) by simplified equations, that disregard the ^-dependence, 
and replace the waveform sech 2 (a/— /ix) by its value at the soliton's center, x = 0: 

l< ^at' = ^ t (1 + 2a) X±1 + 2ax *^ ■ (14) 

An elementary consideration shows that the zero solution of linear equations (jl4j) is double- 
unstable (i.e., through a double eigenvalue) in the region of a < —1/4. 

We checked the stability of the single-component polar soliton (jllj) by means of direct 
simulations of Eqs. (J2J), adding a small random perturbation in the components <f> +1 and 
the value of the perturbation being distributed uniformly between and 0.03 (the same 
random perturbation will be used in simulations of the stability of other solitons, see below). 



The result is that the soliton (fTT|) is always unstable indeed, as an example in Fig. 1. In 
particular, panel (d) in this figure displays the time evolution of 

/+oo p+oo 
\4> {x, t)\ 2 dx, N±i = / |0 ± (x,t)| dx. (15) 
-00 J — oo 

C. Two-component polar soliton 

In the same case as considered above, v — +1, a two-component polar soliton can be 
easily found too: 

O = 0, +1 = ±4>_ l = v / -/isech (y/^JIx) e~* M * (16) 

(note that it again does depend on the value of the nonlinearity coefficient a). In this 
case again, instability is possible due to the parametric gain in linearized equations for small 
perturbations. Indeed, the corresponding equation for the perturbation in the O component 
decouples from the other equations and takes the form 

(id t + fi) Xo = -d^Xo + 2/i sech 2 (sf^x) [(1 + a)x Q ± ax* ] , (17) 

where the perturbation was redefined the same way as in Eq. (fT3|) . O = exp (— ifit) Xo- 
Similar to the case of Eq. (fT2|) . the origin of the instability may be understood, replacing 
Eq. (fTTj) by its simplified version that disregards the x-dependence and substitutes the 
waveform sech 2 (y/—fix) by its value at the central point of the soliton, x = 0: 

z^ = /i[(l + 2a)x ±2oxS], (18) 

cf. Eqs. (I14J) . An elementary consideration demonstrates that the zero solution of Eq. (|TH|) 
(for either sign ±) is unstable in exactly the same region as in the case of Eqs. (jl4j) . a < —1/4. 
Note that the simplified equation (|TH|) illustrates only the qualitative mechanism of the 
parametric instability, while the actual stability border may be different from a = —1/4. 
Indeed, the stability of the soliton (fTBj) was tested by direct simulations of Eqs. (0) with a 
small uniformly distributed random perturbation added to the O component. The result 
shows that the soliton is indeed unstable in the region of \a\ > 1, as shown in Fig. 2, and 
it is stable if \a\ < 1 (not shown here). 

The change of the stability of this soliton at a = +1 can be explained analytically. Indeed, 
splitting the perturbation into real and imaginary parts, Xo( x >^) = Xi( x ,t) + an d 
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looking for a perturbation eigenmode as Xi2( x ^) = e<Tt Ux t 2{x) with the instability growth 
rate (eigenvalue) a, one arrives at a system of real equations [take plus sign in (|I7j)]. 

-aU 2 = -U" - fiUi + 2//(l + 2a)sech 2 (y/^x) U h 
aU x = -U% - //C/ 2 + 2/i sech 2 (V 3 /^) C/ 2 , (19) 

with prime standing for d/dx. The simplest possibility for the onset of instability of the 
soliton is the passage of a from negative values, corresponding to neutral stability, to unstable 
positive values. At the critical (zero-crossing) point, a = 0, equations (fTU|) decouple, and 
each one becomes explicitly solvable, as is commonly known from quantum mechanics. The 
exact solution of the decoupled equations shows that a crosses zero at points a = a n = 
n(n + 3)/4, n = 0,1,2,.... The vanishing of a at the point ao = corresponds not to 
destabilization of the soliton, but just to the fact that Eqs. (pHJ|) become symmetric at this 
point, while the zero crossings at other points indeed imply stability changes. In particular, 
the critical point a\ — 1 explains the destabilization of soliton at a = 1, as observed in 
the simulations, the corresponding eigenfunctions being U\ = sinh (^—fix) ■ sech 2 (y/—fxx), 
U2 = sech (\/— /ix). Critical points corresponding to n > 1 imply, most probably, additional 
destabilizations (the emergence of new unstable eigenmodes) of the already unstable soliton. 
On the other hand, the destabilization of the soliton at a = —1, also observed in the 
simulations, may be explained by a more involved bifurcation, when a pair of the eigenvalues 
a change from purely imaginary to complex ones, developing an unstable (positive) real part 
(a generic route to the instability of the latter type is accounted for by the Hamiltonian Hopf 
bifurcation, which involves a collision between two pairs of imaginary eigenvalues, and their 
transformation into an unstable quartet 



D. Three-component polar solitons 



In the case of v = +1, three-component solitons of the polar type can be found too. One 
of them is 



O = VI - e 2 \/-A* sech (y/^jj,x) e %p,t , 
0+i = -0_i = iev 73 /^ sech (y/^Jix) e"^*, (20) 
where e is an arbitrary parameter taking values —1 < e < +1 (the presence of this parameter 



resembles the feature typical to solitons in the Manakov's system 



23l l26j|). and, as well as 



the one- and two-component polar solitons displayed above, the solution does not explicitly 
depend on the parameter a. We stress that the phase difference of 7r between the <f) +l and 
<p_ 1 components is a necessary ingredient of the solution. 

There is another three-component polar solution similar to the above one (i.e., contain- 
ing the arbitrary parameter e, and independent of a), but with equal phases of the ±1 
components and a phase shift of tt/2 in the O component. This solution is 

O = iVl - e 2 ^/^ sech {^f^jix) e~* M *, 
(f> +1 = 0_ x = iev/ 11 /! sech (^/^jix) e -vrf , (21) 

where the sign ± is the same for both components. 

In addition, there is a species of three-component polar solitons that explicitly depend 
on a: 

O = 1//4 sech (V — f^ x ) 

±1 = x / = £fci sech (y/=fix) e-W, (22) 

where /z ±1 are two arbitrary negative parameters, and the chemical potential is fi = 
— [y + a) (y/—fJ^i + y/—Jj~i) 2 /2, which implies that v + a > (note that this solution 
admits v = —1, i.e., repulsive spin-independent interaction). Each species of the three- 
component polar soliton depends on two arbitrary parameters: either /1 and e [solutions 
(I2TH) and (J2U)], or n_ x and [solution (1221 ]. 

The stability of all the species of these solitons was tested in direct simulations. In the 
case of solutions (|2(J|) and (J21j) . it was enough to seed a small random perturbation only in 
the O component, to observe that both these types are unstable, as shown in Fig. 3 and 4. 

On the contrary, the three-component polar soliton ()22|) is completely stable. This point 
was checked in detail, by seeding small random perturbations in all the components. The 
result, illustrated by Fig. 5, is that the perturbation induces only small oscillations of 
amplitudes of each component of the soliton, and remain robust as long as the simulations 
were run. 

E. Global stability of the solitons 

It follows from the above results that the model demonstrates remarkable multistability, 
as the FM soliton ((7j), the two-component polar solitons (fTIjj) (in the regions of — 1 < a < +1), 
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and the three-component soliton (|22|) may all be stable in one and the same parameter 
region. Then, to identify which solitons are "more stable" and "less stable" , one may fix the 
norm (@J) and compare respective values of the Hamiltonian (jHJ for these solutions, as the 
ground-state solution should correspond to a minimum of H at given N. 

The substitution of the solutions in Eq. (jSJ) reveals another remarkable fact: all the 
solutions which do not explicitly depend on a, i.e., the one-, two-, and three-component 
polar solitons (fTTj) . (JTfiJ . and (|2"0|). (J3TJ, produce identical relations between the chemical 
potential (/i), norm, and Hamiltonian: 



As mentioned above, these relations are different for the single-component FM soliton J7J, 
see Eqs. (jHJ). Finally, for the stable three-component polar soliton (|22jl . the relations between 
fi, N, and H take exactly the same form as Eqs. (jHJ for the FM soliton. 

The comparison of the expressions (|2*3*J and (JHJ leads to a simple conclusion: the single- 
component FM soliton (JJJ and the stable three-component polar soliton (}2*2*|) simultaneously 
provide for the minimum of energy in the case of v = +1 and a > 0, i.e., when both the spin- 
independent and spin-exchange interactions between atoms are attractive (recall a necessary 
condition for the existence of both these types of the solitons is v + a > 0). It is possible 
to state that, in these cases, each of these two species, (J2J and (|2*2J . plays the role of the 
ground state in its own class of the solitons (ferromagnetic and polar ones, respectively). 
Simultaneously, the two-component soliton (JT5J is also stable in the region of < a < 1, as 
shown above, but it corresponds to higher energy, hence it represents a metastable state in 
this region. 

The FM soliton (j2J and three-component polar soliton (|22jl also exist in the case of 
v — — 1 and a > 1, when v + a is positive. In this case too, these two soliton species provide 
for the energy minimum, simply because the other solitons, (JHJ, (fTtjJ. (|2~UJ . and (|2*TJ do not 
exist at z/ = —1. 

For z/ = +1 and —1 < a < (attractive spin-independent and repulsive spin-exchange 
interactions), the Hamiltonian value (J23|) is smaller than the competing one (jHJ- Actually, 
this means that the two-component polar soliton (|TT)j) plays the role of the ground state 
in this only it, among all the polar solitons whose Hamiltonian is given by the 

expression (J2HJ, is dynamically stable (for — 1 < a < +1, see above). As the FM soliton 




(23) 



10 



(J7J) and three-component polar soliton ()22j) also exist and are stable in this region, but 
correspond to greater energy, they are metastable states here. Finally, for v — +1 and 
a < — 1, there are no stable solitons (among the ones considered above). 

The dependence iV(/i) for each solution family provides additional concerning the soliton 
stability. Indeed, the known Vakhitov-Kolokolov (VK) criterion states that a necessary 
stability condition for the soliton family is dM/dfi < 0. It guarantees that the soliton 
cannot be unstable against perturbations with real eigenvalues, but does not say anything 
about oscillatory perturbation modes appertaining to complex eigenvalues. Obviously, both 
relations, given by Eqs. (jHJ) and (J2B1), satisfy the VK condition, i.e., any unstable soliton 
may be unstable only against perturbations that grow in time with oscillations. Indeed, 
numerical results which illustrate the instability of those solitons which are not stable, see 
Figs. 1-4, clearly suggest that the instability, if any, is oscillatory. 



F. Finite-background solitons 



In special cases, it is possible to find exact solutions for solitons sitting on a nonzero 
background. Namely, for v = 1 and a = — 1/2, one can find a two-component polar soliton 
with a continuous- wave (CW) background attached to it, in the following form: O = 0, and 

' 1 

■ i 

For v = a — 1, a three-component polar solution with the background can be found too: 

1 



-ifit 



-ifit 



± sech (v 7- /ix) 
=F sech {\/— /ix) 



(24) 



_ 1 
" 2 

-fie 



-i/it 



V2 



± sech {-\/— fix) 



-ifj,t 



V2 



=F sech (y/— fix) 



(25) 



(in the latter case, the availability of the exact solution is not surprising, as the case of 
v — a — 1 is the exactly integrable one j^] ). 

The stability of solutions (|2*4^) and (|23j) was tested, simultaneously perturbing the O and 
±1 components by small uniformly distributed random perturbations. Figures 6 and 7 show 
that both solutions are unstable, although the character of the instability is different: in 
the former case, the soliton itself (its core) seems stable, while the background appears to 
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be modulationally unstable. In the latter case (Fig. 7), the background is modulationally 
stable, but the core of the soliton (|25|) clearly features an oscillatory instability. 

The modulational instability (MI) of the background in solution (J24)) is obvious, as, even 
without exciting the O field, i.e., within the framework of the reduced equations (jSJ), with 
v + a = 1/2 and v — a = 3/2, any CW solution is obviously subject to MI. As concerns the 
excitation of the O field, an elementary consideration shows that the CW part of solution 
(124)1 exactly corresponds to the threshold of the instability possible through the action of the 
parametric-drive term, 2a0 +1 0_ 1 0Q, in the last equation (J2J). Actually, Fig. 6 demonstrate 
that this parametric instability sets in, which may be explained by conjecturing that the 
initial development of the above-mentioned MI that does not include the <fi Q field drives the 
perturbed system across the threshold of the parametric instability. 



IV. THE DARBOUX TRANSFORMATION AND NONLINEAR DEVELOP- 
MENT OF THE MODULATIONAL INSTABILITY 



In this section, we focus on the integrable case, with v — a — 1, which corresponds to 
the attractive interactions. As explained above, the spinor BEC obeys this condition if a 
special (but physically possible) constraint (JBJ) is imposed on the scattering lengths which 
determine collisions between atoms. Then, Eqs. (J2J) can be rewritten as a 2 x 2 matrix NLS 



equation 



22 



id t Q + d 2 x Q + 2QQ f Q = 0, Q = 



00 <P~l 



(26) 



Equation 



applies 



is a completely integrable system, to which the inverse scattering transform 



2d [i 



[in particular, with regard to the fact that v — a = in the present case, Eqs. (J5J) 
show that system (|2*S|l splits into two uncoupled NLS equations for solutions with <ft = 0]. 
We here aim to develop the Darboux transformation (DT) for Eqs. ()26jl . based on the 
respective Lax pair. 

The Lax-pair representation of Eqs. (|26p is 



U*, * t = V*, 



(27) 
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where U = AJ + P, V = i2X 2 3 + i2XP + iW, with 



. O O Q 

O -I / I -Qt O 



(28) 



^ I- (29) 

Ql -Q f Q 

( *i . 

Here I denotes the 2x2 unit matrix, O is the 2x2 zero matrix, ^ = | is the matrix 

eigenfunction corresponding to A, and ^ 2 being 2x2 matrices, and A is the spectral 
parameter. Accordingly, Eq. (}26|) can be recovered from the compatibility condition of the 
linear system (j27jl . IL - V„. + UV - VU = 0. 

Next, based on the Lax pair ()27j) . and introducing a transformation in the form 

, Ail O . 

* = (A - S)*, S = DAD" 1 , A = , (30) 

O A 2 I 

where D is a nonsingular matrix that satisfies D x = JDA + PD, and letting 

^ = U*, U = AJ + P l5 Pi= ( ° Ql , (31) 

\-Ql o J 

one obtains 

Pi=P + JS-SJ. (32) 
Also, one can verify the following involution property of the above linear equations: if 
^ = ( | is an eigenfunction corresponding to A, where *S/j are 2x2 matrices, then 

is an eigenfunction corresponding to —A*, where we have used the fact that 

Q* = [the matrix Q is symmetric, as per its definition (|2fi)) ] . Thus we can take D in the 
form 

D =i'*'-* ! VA=f AI 

* 2 / V O -A*I 



■ 



to obtain 



S = A + (A + A ) 

O I / Sji -s 22 
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where the matrix elements of S are given by S u = (1 + *i* 2 ^l&Z s i2 = (*2*i 1 + 

vi^;- 1 )- 1 , s 2 i = (^*r x + *i*2 s 22 = (i + ^ 2 ^i i ^i _1 ) 

Finally, the DT for Eq. (}2T)|) follows from relation (}3*2"j) : 



It should be noted that = S12 implies the above-mentioned property, = Q{. From 
Eqs. ()33|) and (|27|). it can be deduced that Eq. generates a new solution Q x for Eq. 
(I26j) once a seed solution Q is known. In particular, a one-soliton solution can be generated 
if the seed is a trivial (zero) state. Next, taking Qx as the new seed solution, one can derive 
from Eq. (}3*3*j) the corresponding two-soliton solution. This can be continued as a recursion 
procedure generating multi-soliton solutions. 

In the following, we will be dealing with solutions to Eq. (}2l)j) under nonvanishing bound- 
ary conditions. The simplest among them is the CW solution, with constant densities of the 
three spin components, 



where a c and (3 C are real constants, and A; is a wave number. In this solution, the constant 
densities are equal in components ±1 , while the phase difference between the components 

is 7T. 

Applying the DT to the CW solution Q c , solving Eqs. (|2*7|) for this case, and employing 
relation (}3*3*j) . we obtain a new family of solutions of Eq. (|2l)Jl in the form 



Qx = Q + 2(A + A*)S 



(33) 




<p c = kx + [2(a 2 c + Pi) - k 2 ] t, 



(34) 



Q 1 = [A c + 4^(I + AA*)- 1 A]e^, 



(35) 



where the following definitions are used: 



A = (Ue e -^ + «r 1 A c )(K- 1 A c IIe < '- i<p + I) 



-1 



(36) 



<p = M R x - [2fMj + (k + 2r])M R }t, 



9 = Mix + {2£M R - (k + 2r t )M I }t, 



(37) 
(38) 




(39) 



14 



.(3a. 

~(zk — 2X + iM), X = £ + iri is the spectral parameter, and II = | | is an arbitrary 

a 7 

complex symmetric matrix. Solution ()35|) reduces back to the CW (|34|) when £ = 0. It is 
worth noting that the three-component polar soliton (J25|) considered in the previous section 
is not a special example of the solution (JSHJ), because the background densities in solution 
(P5|) have equal phases in components ±1 . 

In particular, with A c = O (zero background), Eq. (|33j) yields a soliton solution, that 
can be written as 

4£[IIie- es + (cr^cr^e 93 det ni]e^ 
Ql ~ e- 2 ^ + l + e 2e Hdetn 1 | 2 ' ^ ' 

where 9 S = 2£(x — Ar/t) — 6 , <p a = 2r/x + 4(£ 2 — r/ 2 )t, and 8 is an arbitrary real constant 

which determines the initial position of the soliton, cr y is the Pauli matrix, and Hi is the 

polarization matrix j^ . 

n 1= (2M 2 + |/?| 2 + | 7 | 2 )- 1/2 n 

Solitons (|4(jp are tantamount to ones found in Ref. |22j, where they were classified into 
the above-mentioned distinct types, the FM and polar ones. Indeed, one can see from 
expression (|4T)j) that, when detlli = 0, the soliton (|4*0|l represents a FM state in the spinor 
model. However, when det Hi 0, the solution ()40|) may have two peaks, corresponding to 
a polar state (the smaller detlli, the larger the distance between the peaks). 

We will now produce another particular example of exact solutions (|35jl . that describes 
the onset and nonlinear development of the MI (modulational instability) of CW states. 
Indeed, under the conditions k = 2r] and a 2 + (3 2 c > £ 2 , Eq. ()39|) yields M\ = 0, hence 9 
is independent of x [see Eq. (J37J)], and the solution features no localization. Instead, it is 
x-periodic, through the x-dependence in ip, as given by Eq. in this connection, we 

notice that 

Ml = 4 (« 2 + (3 2 c ) - Ai 2 (41) 

does not vanish when Mj = 0. This spatially periodic state may actually display a particular 
mode of the nonlinear development of the MI, in an exact form. To analyze this issue in 
detail, we consider a special case with II = eE, where E is a matrix with all elements equal 
to 1, and e is a small amplitude of the initial perturbation added to a CW to initiate the 
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onset of the MI. Indeed, the linearization of the initial profile of solution ()35|) in e yields 



where p = 1 - £(2£ + iM R )/{a 2 c + (? c ) with \p\ = 1, Xl = £M*(2^ - M fl )/(a 2 + /3 2 ), 



X 2 = /? c Xi/( a c + /^c)> Mr = 2\/ a 2 + /? 2 — £ 2 , and cr 2,:r are the Pauli matrices. Clearly, 
the first term in Eq. ()42j) represents the CW, while the second and third ones are small 
spatially-modulated perturbations. 

Comparing the exact solution (jHSJ) obtained in this special case with results of direct 
numerical simulations of Eqs. (j2J) with the initial condition (}4*2"j) (not shown here), one can 
verify that the numerical solution is very close to the analytical one, both displaying the 
development of the MI initiated by the small modulational perturbation in the waveform 
(jl2| - Figure 8 shows the result in terms of the densities of the three components, as provided 
by the exact solution (|35|). From this picture we conclude that atoms are periodically 
transferred, in the spinor BEC, from the spin state mp = into ones mp = ±1, and vice 
versa. 

As said above, the exact solution (}35|) to Eqs. (j2j) is only valid under the special condition 
a = 1 (i.e., 2g = — (fe)- In the general case (a ^ 1), the onset of the MI can be analyzed 
directly from equations (J2J) linearized for small perturbations, as it was actually done in Ref. 
[is| . We will briefly recapitulate the analysis here (in particular, to compare results with 
the above exact solution). To this end, we take the CW solution (J34J) (with k = 0, which 
can be fixed by means of the Galilean transformation), and consider its perturbed version, 



where b±i and b are weak perturbations. Substituting this expression into Eqs. (0), lin- 
earizing them, and assuming the perturbations in the usual form, 



where frj-Rj are real amplitudes, the linearized equations give rise to two branches of the 
dispersion relation between the wavenumber K and frequency uj of the perturbation, 



Qi(x, 0) w [A c p - eXiEcos(M R x) - ex 2 cr z cr x A c e iMRX ] exp(ikx) 



(42) 




bj(x, t) = bj R cos(Kx — cut) + ibji sin(i^x — cut), 



j = +1,0,-1, 



u? = K 2 (K 2 - Aao? c - 4a/? 2 ); 



(43) 
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u 2 = K (K - 4va 2 c ~ We) • ( 44 ) 

The MI sets in when uj 2 given by either expression ceases to be real and positive. Obvi- 
ously, this happens if, at least, one condition 

K 2 < Aa(a 2 c + p 2 c ), K 2 < 4u(a 2 c + (3 2 c ) (45) 

holds. Note that these conditions agree with the above exact result, which corresponds to 
a — v — 1. Indeed, in that case we may identify K = Mr, according to Eq. (j42j) . and then 
Eq. ()41|) entails the constraint, < 4(a^ + /3^), which is identical to conditions ()45|) . 
Note that the MI may occur in the case of the repulsive spin-independent interactions in 
Eqs. (J2J), i.e., if v — — 1, provided that the nonlinear-coupling constant a, which accounts 
for the strength of the spin-dependent interaction, is positive. 

We also explored stability of the soliton solutions (|4T)|) . of both the polar and FM types, 
against finite random variations of the coupling constant a in time. The issue is relevant, 
as the exact solutions are only available for a = 1, hence it is necessary to understand 
if the solitons survive random deviations of a from this special value corresponding to the 
integrability (in fact, this way one may test the structural stability of the solitons generated in 
the exact form by the inverse-scattering/Darboux transform). Physically, this situation may 
correspond to a case when the scattering length, that determines the nonlinear coefficient a 
in Eqs. (0), follows random variations of an external magnetic Q or laser field which 
affects the scattering length through the Feshbach resonance. The evolution of the density 
profiles in the thus perturbed solitons is displayed in Fig. 9. This figure and many other 
runs of the simulations suggest that the solitons of both the FM and polar types are robust 
against random changes of a, once they are dynamically stable as exact solutions to the 
integrable version of the model. Note also that the FM soliton seems more robust than its 
polar counterpart. 



V. CONCLUSION 



In this work, we have presented new soliton states in the model based on the coupled NLS 
equations describing the dynamics of F = 1 spinor BECs. In the general (non-integrable) 
version of the model, several types of elementary exact soliton solutions were reported, which 
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include one, two, or three components, and represent polar or ferromagnetic (FM) solitons. 
Their stability was checked by direct simulations, and, in some cases, in an exact analytical 
form, based on linearized equations for small perturbations. The global stability of the 
solitons was analyzed by comparing the respective values of the energy for a fixed norm 
(number of atoms). As a result, we have found stable ground-state solitons of the polar and 
FM types; stable metastable solitons coexisting with the ground-state ones are possible too, 
in certain parameter regions. 

In the special case of the integrable spinor-BEC model, we have applied the Darboux 
transformation (DT) to derive a family of exact solutions on top of the CW state. This 
family includes a solution that explicitly displays full nonlinear evolution of the modulational 
instability (MI) of the CW state, seeded by a small spatially periodic perturbation. Finally, 
by means of direct simulations, we have established robustness of one-soliton solutions, of 
both the FM and polar types, under finite perturbations of the coupling constant. In fact, 
the latter result implies the structural stability of these solitons. 

This work was supported by the NSF of China under grants 60490280, 90403034, 
90406017, and Overseas Scholar Foundation of Shanxi Province. The work of B.A.M. was 
supported, in part, by the Israel Science Foundation through a grant No. 8006/03. 
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Figure Captions 

Figure 1. Evolution of the single-component polar soliton (jllj) with a small random 
perturbation added, at t — 0, to the <p +l and <ft_ l components. Parameters are v — 1, 
a = —0.5, and fi = — 1. 

Figure 2. The solid, dashed, and dotted curves show, respectively, the evolution of the 
norms N , N± of the O and ±1 components, and total norm iV [see Eqs. (II 5|) and (@J)] in 
the two-component polar soliton ()16|) perturbed by a small random perturbation introduced 
in the O component. Parameters are v — 1, \i — —1, and (a) a = 1.5 for <p +1 (f)*L l < 0; (b) 
a = 1.5 for > 0; (c) a = —1.5 for 4> +1 (p*_ 1 < 0; (d) a = —1.5 for > 0. 

Figure 3. The evolution of the three- component polar soliton (j2U|) under the action of 
a small random perturbation initially added to the O component. Parameter are v — 1, 
fi = —0.8, e = 2/y / 5, and (a) a = 0.5; (b) a = 1.5; (c) a = —0.5; (d) a = —1.5. The meaning 
of the solid, dashed, and dotted curves is the same as in Fig. 2. 

Figure 4. The same as in Fig. 3 for the three-component polar soliton (}2"Tj) . Parameters 
are u — 1, /j, — —0.8, e = 2/v^5, and (a) a = 0.5; (b) a = 1.5; (c) a = —0.5; (d) a = —1.5. 
The meaning of the solid, dashed, and dotted curves is the same as in Figs. 2 and 3. 

Figure 5. The evolution of the three-component polar soliton f!22j) with initially added 
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random perturbations. Parameters are v = 1, = —1.21, = —0.25, and (a,c) 
a = 0.5, (b,d) a = —0.5. In cases (a) and (b), the random perturbation was added to the 
0g component, and in cases (c) and (d) - to components ±1 . The curves marked by 1, 
2, 3, and 4 depict, respectively, the evolution of the norms iV +1 , iV_i, iVo and N, see the 
definitions in Eqs. (fT5|) and (jJJ). 

Figure 6. The evolution of the two-component polar soliton on the CW background, 
(J2H), under the action of initial random perturbations added to the O and ±1 components. 
Parameters are v = +1, a = —1/2, and \i = — 1. 

Figure 7. The same as in Fig. 6, for the finite-background soliton (|25j) . Parameters are 
z/ = a = 1, and fi = —0.36. 

Figure 8. The nonlinear development of the modulation instability, as per the exact 
solution (j33j) under the conditions k = 2r] and a 2 + (3 2 c > £ 2 . Parameters are k = 0.6, 
r] = 0.3, £ = 1, a = P = j = e~ 8 , a c = 1.2, and /? c = 0. 

Figure 9. The evolution of the density distribution in the soliton obtained from solution 
(f33j) by adding a randomly time-dependent perturbation of the nonlinear-coupling constant 
a, with the perturbation amplitude ±5%. Parameters are v = a = 1 (as concerns the 
unperturbed value of a), £ = 1, rj = —0.03, and 6 = 0. (a) The FM soliton, with cti = 0.48, 
(5 X and 7 X being determined by conditions a\ = /3 1 7 1 and 2 |«i| 2 + \ f3 l | 2 + I7J 2 = 1; (b) the 
polar soliton, with a\ = 0.53 and / 9 1 = 0.43, r y l being determined by 2 |ai| 2 + |/3i| 2 + |7i| 2 = 1. 
The solid and dotted curves show, respectively, the exact solutions and the perturbed ones 
produced by numerical simulations. 
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